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Abstract. We present a Lattice-Boltzmann method for simulating self-propelling ( active) colloidal particles 
in two-dimensions. Active particles with symmetric and asymmetric force distribution on its surface are 
considered. The velocity field generated by a single active particle, changing its orientation randomly, and 
the different time scales involved are characterized in detail. The steady state speed distribution in the 
fluid, resulting from the activity, is shown to deviate considerably from the equilibrium distribution. 
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1 Introduction 

Motion at low Reynolds number is an intriguing phenom- 
ena. How does a micron sized bacteria moving in a fluid 
change the velocity field of the fluid around it? Does this 
local stirring change the viscocity and temperature of the 
fluid? Such questions have gained importance, thanks to 
, recent experiments on the motion of bacteria in two di- 
mensional fluid films Theoretical work on active 
particle suspensions has addressed the organization of ac- 

■ tive particles and their collective motion [HE], as well as 
pattern formation induced by activity . The peculiar- 
ities of the rheological response of active suspension has 
also been analyzed [7||H]. 

The collective, hydrodynamic behavior of active inter- 

■ acting organisms has been analyzed disregarding the dy- 
namics of the embedding solvent U- In doing so, the mi- 

, croscopic details of propulsion and particle interactions 
were hidden in the phcnomenological coefficients. In a 
more microscopic description, Ramaswamy and co work- 
ers model the active particles as point force dipoles [H] 
and study the possibility of nematic ordering in such a 
system. More complicated force distributions on the ac- 
tive particles are difficult to carry out in theoretical treat- 
ments. In this paper we present a Lattice-Boltzmann sim- 
ulation scheme for active colloidal particles. The model 
can be used to study the rheological behaviour and self- 
organization in active particle suspensions 9 . 

The Lattice-Boltzmann (LB) method relies on the fact 
that the macroscopic dynamics of fluid flow does not de- 
pend on the microscopic details of the solvent. LB is a 
mesoscopic kinetic model that captures the collective macro- 



scopic dynamics of fluid flow. The method involves the 
discretization of both space and time, thereby reducing 
the number of degrees of freedom of the system, which 
in turn leads to a reduced computational effort. By intro- 
ducing simplified kinetic rules one captures the essential 
macroscopic behaviour to a high degree of accuracy. 

Over the years, LB has been used to simulate a wide 
variety of systems It provides a faithful discretization 
of Navier-Stokes equation near incompressible-fluid flow, 
and has been used for large scale fluid dynamic simula- 
tions m. It has also proved to be a flexible method for 
the study of complex fluid flows and has been applied to 
study porous medium flows, binary fluid separation kinet- 
ics as well as colloidal and polymer suspensions ^2] ■ 

Among the different proposals to simulate rigid parti- 
cles suspended in a fluid (such as a colloidal suspension) 
in LB we will follow the method introduced by Ladd 
which defines the solid through the links that join particle 
and fluid nodes, called the boundary links. The suspended 
objects interact with the surrounding fluid through the 
conditions specified at the boundary links. In the original 
version, the fluid also occupies the interior of the closed 
surface defined by the boundary links. The inertia of the 
inner fluid may affect the dynamics of suspended colloids 
at short time scales. More recently, Nguyen and Ladd |14j . 
have introduced a variant of this method which excludes 
the fluid at the interior of the colloid. While for single 
phase fluids such a variant is not crucial, excluding the 
fluid from inside the colloids cannot be avoided in the case 
of fluid mixtures [15j . where the concentration gradients 
in the neighbourhood of the colloid must be accounted for 
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appropriately. In both descriptions the appropriate stick 
boundary conditions of the fluid at the solid surface are 
accomplished by bouncing back the incoming fluid densi- 
ties in the frame of reference of the moving particle. 

In this paper we present the first two-dimensional LB 
model of a suspension in which the solid particles are 
driven by internal forces. In section II, we describe in de- 
tail the method used to simulate the active motion of the 
particles through appropriate multipolar force distribu- 
tions. Section III concentrates on the results for a sin- 
gle active particle, in the absence of thermal fluctuations, 
when the force distribution does not induce any net mo- 
tion - a " shaker" . Section IV describes the peculiarities of 
a particle that can displace - a " mover" . We conclude with 
section V, where we discuss the main results obtained and 
indicate future directions. 



2 Active Colloid Model 

The LB approach has been widely used to study the dy- 
namics of colloidal suspensions ^2] ■ We use the algorithm 
proposed by Ladd for a colloidal suspension and extend it 
to incorporate active particles. In the following section we 
briefly describe the discretized LB approach. 

2.1 Fluid Lattice Model 

The central dynamic quantity in LB is the one particle 
distribution function, 7ij(r, Cj, t), which corresponds to the 
density of fluid particles at the lattice site r moving with 
velocity c, . All the length and time scales used are in terms 
of lattice units, which implies that the particles move a 
unit distance per unit time. 

Having defined the distribution function, rij(r, Cj,t), 
one can recover the hydrodynamic fields as its various mo- 
ments. In particular, the zcroth moment gives the local 
density of the fluid p = m, the local momentum of the 
fluid j = pu = Y). JijCj. while the momentum flux is given 
by the second moment, II = Yli n i c i c ii where CjCj is a 
diadic product. 

In LB, at each time step the distribution at every lat- 
tice site is updated in two substcps; a propagation and a 
collision event. In the propagation event, is transferred 
to the corresponding neighbouring site along the link cor- 
responding to its velocity, Cj. In the collision substep, the 
distribution at every site is relaxed toward its local equi- 
librium distribution in a way that ensures local mass and 
momentum conservation. The overall dynamics can be ex- 
pressed as 

n ?; (r + c l; t+ 1) = n;(r,i) -nf ), (1) 

3 

where the first term in the r.h.s. expresses the propa- 
gation, and the second term corresponds to collisions. £ij 
is the linearized collision matrix while rij 9 is the equilib- 
rium distribution which depends on the local values of p 



and u. For simplicity's sake, we restrict ourselves to the 
particular case in which 77 relaxes to its equilibrium value 
in one time step. In this case the equilibrium distribution 
function can be simplified to 

nJ 9 = p[ag'+o?u-c i ] ) (2) 

with fluid viscosity set to r\ = The coefficients in 
the equilibrium distribution function are obtained impos- 
ing isotropy, mass and momentum conservation. The post- 
collision distribution can now be expressed as 

n 4 + A l {n l ) = a^'p + a c {j a c la (3) 

for zero Reynolds number flows. The results obtained in 
this paper correspond to a D2Q7 lattice, which means 
a two dimensional triangular lattice, with seven velocity 
directions at every site, see figure ^ For this lattice, the 
amplitudes of the equilibrium distribution are 

a° = l -2^,^ = ^/3,^ = 1/3, (4) 

where c s is the speed of sound in the fluid. Since the dis- 
tribution function must be positive, c s is restricted to 
< c s < -7=. For convenience, the density of particles 
at every point moving in each direction is set to unity, 
leading to p eq = 7 at every lattice site. 




Fig. 1. The seven velocity directions a, where i — 0, 1...6. 



2.2 Solid particles 

We follow Ladd's approach and define the solid particles 
through the set of boundary links which define a closed 
surface. These boundary links join neighbouring fluid and 
solid nodes. This procedure allows for simulation objects 
of any shape, as shown in figurc|21 As it moves, the center 
of mass of the particle makes a smooth trajectory in space. 
It is worth noting that the shape of the object changes 
as it moves, since the number of boundary links will in 
general change. Such lattice artifacts can be systematically 
reduced by increasing the particle's size. 

Apart from the propagation and relaxation events which 
define the LB dynamics, the density of fluid particles mov- 
ing along a boundary link are bounced-back to enforce 
stick boundary conditions as described below. Lattice nodes 
on either side of the boundary surface are treated on the 
same footing, so that the fluid fills the space inside and 
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which give rise to net force and torques acting on the 
particle, 

P(r,t+1) = ^2 [m(T,t) -n v (T + a,t) 

-2oiCi/o(r)ui, • Ci] Cj (10) 

T(v b ,t+l) = Y,nxF b (r b ,t+l) (11) 



Fig. 2. Discrete representation of a solid particle along with 
the boundary links, when its center of mass is displaced from 
the lattice node. The dashed circle represents the particle. 



where the sums run over all the boundary links. Then the 
particle linear velocity and angular velocity are updated 



using 



U(t+1) = U{t) + 2M- 1 F(t + l) (12) 



outside the solid particles; the interior fluid being kept 
only for computational convenience. If needed, there are 
alternative approaches which exclude the interior fluid ^3] . 

At each boundary link there are two distributions rii (r, t) 
and rii' (r + Cj, t), corresponding to velocities Cj and Cj/ 
(cj/ = — Cj) parallel to the boundary link connecting r, 
and Ti+Ci. If the boundary is stationary, the incoming dis- 
tributions are simply reflected back in the direction they 
come from. If the boundary is moving, upon reflection, a 
part of the momentum is transferred across the boundary 
node. These boundary node update rules can be written 
as 

n,(r + c 4 ,i) = m>(r + Ci,t) + 2a^pu b • c, (5) 



n(t + i) = n{t) + 2r 1 r{t + i) (13) 

where M is the mass of the particle and I is its moment 
of inertia. 



rii (r j/ , t) = rii (r , t) — 2al* pu b ■ a 



(0) 



where the boundary link velocity, u b , is calculated by 



u b = U + n x (r + y -R). 



(7) 



Here U is the center of mass velocity of the particle, f2 is 
the angular velocity, r + Cj/2 = r b is the position of the 
boundary link while R is center of mass position. 

Due to the bounce back, the fluid exerts a net force 
on the particle. As a result, there is an exchange of mo- 
mentum at the boundary, but the total momentum of the 
particle and the fluid combined is conserved. Using this 
condition, the net force and momentum acting on the par- 
ticle can be calculated. The force and the torque at every 
node can be expressed as, 



F b (r b ) = 2 K(r,i) - ?v (r + c t , i) 

-2aip(r)u 6 • Cj] c t (8) 
r b = v b x ~F b {r b ) (9) 



2.3 Self-propulsion 



So far we have been discussing colloidal particles that are 
passive. In presence of thermal fluctuations the boundary 
interaction described above will ensure that the particles 
execute Brownian motion. We will now modify the above 
algorithm in order to incorporate the self-propulsion of 
the particles, their "activity". Since by Newton's third 
law, the net force acting on a self-propelled particle must 
be zero, we thus have to distribute forces, at the bound- 
ary nodes, such that there is no net force exerted by the 
particle. To this end, we divide the particle into two parts 
and distribute equal and opposite forces in the two halves. 
The simplest form for such a self-propulsion is a force 
dipole |T7j. 

To implement this scheme, first a random direction 
is chosen, denoted by the unit vector h. The particle is 
then divided, perpendicular to this axis, into two parts 
at a distance d from its geometric center. We then assign 
forces to each of the boundary nodes such that the net 
force in one sub-domain is equal and opposite to that in 
the other sub-domain. Let us consider a more specific case 
where the net force exerted by the particle through one 



sub-domain is F , given by 



F n 



(14) 



Since the particle can only interact with the fluid through 
the boundary links, it can exert force only along those di- 
rections. This force along the boundary links should be 
chosen such that their sum is equal to F a . Obviously, 
there are many ways of distributing this force. We chose a 
scheme where the force fi along each of these links i have 
components given by 



fia 



(15) 
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where the sum is over all boundary links in the sub-domain 
that has a positive component along F a . 

The resulting force distribution is shown, schemati- 
cally in figure These are the elementary force that will 
enter into the LB method. 




Fig. 3. Schematic representation of the force distribution on 
the particle. The dashed line is the plane which divide the 
particle into two. The sum of the forces in each part is equal 
and opposite. The filled rectangle is the geometric center and 
the filled circles are the force centers in each half. 



The force acting on each boundary node is then 

I','. = fi*e x + f iy e y (16) 

This force will result in an additional contribution to 
the boundary node velocity which can be obtained by in- 
verting canllOl 



2 a\p 



m(r,t+l) 



'(r + Ci,t + l) 



2c 



(17) 

The boundary node velocity that will enter into the 
calculation of force in eqn [3] and torque in eqn [U] is now 
the sum of active and passive contributions, + u^. The 
procedure is summarized in figure 0] 

We restrict ourselves to the zero temperature case, 
wherein the fluid velocity field is generated solely by the 
force exerted by the active particles on the fluid. The effect 
of the parameter d, the distance of the partitioning line 
from the geometric center of the particle, is to change the 
number of boundary nodes along which the internal forces 
are acting. When d ^ the number of terms in the sum 
in eqn^]will be different for the two parts of the particle. 
We will now discuss the implication of this asymmetry in 
force distribution. 



2.4 Movers and Shakers 

As we discussed before since the forces exerted by the 
colloid are internal, the net force should be zero. This 
means that the forces exerted by the boundary nodes on 



Direction h 



Divide perp. to h 









F c + F / q 




ax x a y y 




Fig. 4. Summary of the procedure for distributing forces on 
the particle 



the fluid, to lowest order, can only be a force dipole |17| . 
However, the net force exerted by the fluid on the particle 
depends on the symmetry of this force distribution. 

In our simulations, in each part of the particle, we can 
define the location of a force center as 



R 



fc 



E^ 2 



(18) 



where the summation is over all the boundary nodes in 
that segment. If b and b' are the distances to this force 
center from the geometric center of the particle in each 
half (sec figure |3J)), then w = \b — b'\ is a measure of the 
asymmetry in the force distribution. The quantities w and 
d are then linearly related. 

When w is non zero, the velocity field in the fluid in 
the front and back segments of the particle are different. 
This asymmetry in the velocity field results in a net force 
on the particle which generates its motion. We call this 
a " mover" . For single particles in a symmetric environ- 
ment this asymmetry in the velocity field vanishes when 
w is zero. However the particle still generates a non-zero 
velocity field in the fluid, in other words it behaves as a 
" shaker" . Stated differently, if the division of the particle 
is about the geometric center (d = 0), the force centers 
in each part are at equal distances from the center, and 
one obtains an apolar particle, otherwise it is polar |7j. 
Both polar (mover) and apolar (shaker) particles disturb 
the fluid around it. The apolar particle by symmetry can- 
not move and will only generate a symmetric flow around 
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it (see figure EJ). But the polar particle will get displaced 
as a result of the applied force distribution, since it can 
induce an asymmetric flow field around it (see figure |5J). 



3 Results - Single Shaker 

We will first look at the case of symmetric force distri- 
bution on a particle placed at the center of the cell 1 . In 
this case the center of mass of the active particle does not 
show any motion along any preferential direction. As a 
result of the force applied by it, the particle will gener- 
ate a steady velocity profile around it. In order to sample 
the statistical properties of the velocity field, we randomly 
change the direction along which force is applied at reg- 
ular intervals, we call this the switching period, and the 
time interval between them as the switching time. We have 
found that the maximum fluid disturbance is achieved if 
we switch the direction in the time scale in which the fluid 
flow adapts to the perturbation imposed by the particle, 
which scales as l?jv. We obtain this time scale as the 
time taken for the average fluid velocity, to reach a steady 
value. In figure we plot the time derivative of the nor- 
malized average speed of the fluid particles as a function 
of time for different values of the system size. As can be 
seen from the figure, good data collapse is obtained if we 
scale the time axis by L 2 . For the rest of this paper we 
measure time in units of L 2 jv . 

The velocity profile around the particle at this steady 
state condition is shown in figure The force applied 
by the particle on the fluid creates four vortices that are 
symmetrically placed. As a result, there is no net force 
acting on the particle and it does not get displaced. In fig- 
ure 0we show the fluid velocity distribution for different 
switching rates. It is clear from the figure that there is a 
dependence on the switching rate. Maximum disturbance 
in the fluid is created when the switching rate is roughly 
\/{10L 2 /v) 1 which corresponds to a switching time of the 
order the time needed to reach the steady state . The ve- 
locity distribution in the steady departs clearly from a 
Maxwellian. This is expected given the fact that the sys- 
tem is clearly out of equilibrium. In figuredwe compare a 
Maxwellian with the distribution obtained from the sim- 
ulations to show that except for very slow switching rates 
we have significant deviation from it at large velocities. 



4 Results - Single Mover 

When the force centers are not symmetrically placed about 
the particle's geometric center ( w ^ 0) the interaction of 
the particle with the fluid results in a net force on the par- 
ticle. The particle now has a direction set by the vector 
d, along which it starts to move. Starting from rest, the 

1 We work with two values,Af = 1000 and 10000 for the 
mass of the particle, which corresponds to particle densities, 
p p — 12/0/ and p v = 120p/ , respectively. The simulations are 
stable for both these values of the mass. However for values of 
M < 1000 the simulations are unstable. 




ln(x v/L 2 ) 



Fig. 5. Time derivative of the average speed of the fluid as a 
function of scaled time, in the presence of a single shaker, v* is 
the average speed at steady state. Parameters used are R = 5, 
M = 1000 . All the system sizes simulated fall on the same 
curve, showing a relaxation with two algebraic decay. 



Fig. 6. The fluid velocity profile around a shaker. R=5, L=100, 
and Af=1000 

particle and the fluid reach a steady state velocity on a 
time scale L 2 /v, as was the case for the shaker. However 
the fluid velocity profile, at this steady state, is quite dif- 
ferent from that of the shaker. As can be seen in figure 0D 
the vortices in the front and back of the particle are now 
asymmetric, clearly indicating the direction of motion of 
the particle. 

We have also analyzed the decay of the particle's veloc- 
ity if its activity is switched off. For the systems analyzed, 
such a decay is consistent with a double exponential. This 
is depicted in the figure|!|] The time scales n and r 2 origi- 
nate from different dissipation mechanisms. The first time 
scale comes from the friction between the particle and the 
surrounding fluid, and is the result of particle's inertia. 
In a short time after activity is switched off, the velocity 
profile changes to one with no slip at the boundary. Fur- 
ther dissipation is due to the drag resulting from the dis- 




6 



S. Ramachandran, P.B.Sunil Kumar and I. Pagonabarraga: Lattice-Boltzmann for active colloids 




Fig. 7. Normalized fluid velocity distribution P(v') vs v' for 
different switching times in a semi-log plot for a shaker, v'is 
the velocity scaled by the average speed of the fluid velocity. 
Other parameters are 7?=5, £=100 and M=10,000 
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Fig. 8. Asymmetric velocity profile of the fluid around a mover 
at steady state. R=5, L=W0 and M=1,000 



tortions in the fluid created by the particle motion. These 
distortions travel through the fluid and are eventually dis- 
sipated. In two dimensions, the resistance offered to the 
particle motion due to this Stokes' drag, has a weak depen- 
dence on both the system size and the particle size |18| . 

The velocity of the particle depends on two factors, the 
strength of the force in each half (/) and the asymmetry 
in the force center positions with respect to the geometric 
center (w). The dependence of the particle velocity on w 
is shown in figure ITUl 

As happens with the shaker, the activity of the particle 
results in a non-Maxwcllian distribution for the fluid node 
velocities. This is shown in figure lTTI This distribution also 
has a dependence on the switching time of the particle, 
with the maximum velocity obtained when the switching 
time is comparable the the diffusion time set by L 2 jv. The 
curves are analogous to those obtained in fig. EI except for 
the fact that the tail is longer here. 




Fig. 9. The velocity of the colloidal particle as a function 
of time after the driving is switched off. The continuos line 
is a fit to the function of the form v\e~ l + vie~ l ^ T2 with 
ri = 0.0016 , t 2 = 0.01 , ui = 0.8 and v 2 = 0.2. The data is for 
a particle of radius R = 7, mass M = 1000 . The size of the 
simulation box is L — 100. The velocity has been scaled by the 
particle velocity at the time when the activity was switched off 
i.e.at t' = 




1.5 



& 1 



0.5 







w 

Fig. 10. Speed of the particle for different values of w. Other 
parameters are R=5, L=100 and M=1000. Here (v' p ) is the 
average particle particle velocity in units of Ljv 



4.1 Many Particles 

We will now look at the effect that a suspension of ac- 
tive particles have on the fluid. We restrict ourselves to 
the dilute limit where direct interactions have a negligible 
contribution. 

The fluid velocity distribution with five active parti- 
cles in a 100 x 100 lattice is presented in figure ^| for 
both movers and shakers . The non-Maxwellian nature of 
the distribution is highlighted here by comparing the p(v) 
for large v with that obtained from the best fit Maxwell 
distribution. 

In order to further check the effect of shakers and 
movers on the fluid we compute the diffusion coefficient 
of a passive bead, having the same radius as the active 
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100 



Fig. 11. Normalized velocity distribution P(v ) vs v of the 
fluid nodes for different switching times (shown with arrow 
marks) for a mover. Velocity, v' , is the velocity scaled by the 
average particle velocity. Other parameters are R=5, L—100 
and M=10,000 
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Fig. 12. P(v') vs v' for five movers and shakers in a semi-log 
plot. Other parameters are the same as in 1111 
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Fig. 13. ^-77-^ vs t with five active particles and one pas- 
sive particle. The position of the passive particle is tracked 
as a function of time to obtain the data. R=5, L=100 and 
M=10,000 



results in a shaker particle, which disturbs the fluid around 
it by producing a velocity field but without producing any 
net displacement. A mover is a result of asymmetric force 
distribution. These particles can move by using the asym- 
metry in the velocity field around it. We show that the 
velocity of the mover is directly proportional to the asym- 
metry in the position of the force distribution. 

The statistical properties of the velocity field generated 
by the activity is shown to deviate from that of equilib- 
rium fluids, with the velocity distribution having a gener- 
alized exponential exponential tail p(v) exp(— (v/vo) a ■ It 
is shown that the momentum diffusion and the different 
dissipative mechanisms are consistent with what is well 
known in two dimensional fluids. 

The model could be used to study the emergence of 
collective behavior in self propelled organisms as well as 
the rheological properties of the suspension. These aspects 
are currently under investigation. 



ones, suspended in the fluid along with several other active 
particles. The mean square displacement of this particle 
as a function of time is shown in figure HHl During early 
times, the (r 2 ) shows ballistic behavior, and at later times 
a crossover regime leads eventually to a purely diffusive 
regime. 

The above two measurement indicate that the movers 
keep the system at a higher temperature than the shakers. 



5 Conclusion 

We have presented a model for the numerical simulation of 
self-propelled colloidal suspensions based on the Lattice- 
Boltzmann method for colloids. By devising a way of adding 
internal forces to the boundary interaction between the 
fluid and the solid particle so as to satisfy the net zero 
force condition, we are able to model two types of parti- 
cles - movers and shakers. A symmetric force distribution 
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